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We present a novel framework to decompose three-nucleon forces in a momentum space partial- 
wave basis. The new approach is computationally much more efficient than previous methods and 
opens the way to ab initio studies of few-nucleon scattering processes, nuclei and nuclear matter 
based on higher-order chiral 3N forces. We use the new framework to calculate matrix elements of 
chiral three-nucleon forces at N 2 LO and N 3 LO in large basis spaces and carry out benchmark calcu¬ 
lations for neutron matter and symmetric nuclear matter. We also study the size of the individual 
three-nucleon force contributions for 3 H. For nonlocal regulators, we find that the sub-leading terms, 
which have been neglected in most calculations so far, provide important contributions. All matrix 
elements are calculated and stored in a user-friendly way, such that values of low-energy constants 
as well as the form of regulator functions can be chosen freely. 

PACS numbers: 21.30.-x, 21.45.Ff, 13.75.Cs 


I. INTRODUCTION 

The importance of chiral three-nucleon (3N) forces has 
been demonstrated in numerous microscopic calculations 
of few-body scattering processes, nuclei and nuclear mat¬ 
ter, see mm for recent review articles. Contributions 
to 3N forces (3NFs) start to appear at next-to-next- 
to-leading-order (N 2 LO) in the chiral expansion within 
Weinberg’s power counting framework 011, whereas 
the subleading chiral 3NFs at next-to-next-to-next-to- 
leading-order (N 3 LO) have the particular feature that 
they do not include any new low-energy constants 0 - 
0. So far, all microscopic calculations for finite nuclei of 
mass A > 3 based on chiral EFT interactions were limited 
to 3NFs at leading order (N 2 LO), whereas NN forces are 
commonly included up to N 3 LO, see m for a fifth order 
analysis of NN scattering. For NN forces, it is known that 
N 3 LO contributions are important for a high-precision fit 
of scattering phase shifts and mixing angles up to labo¬ 
ratory energies of E ~ 200 MeV 0CLM5]. It is clearly 
desirable to extend these studies and perform complete 
N 3 LO calculations, which would require the inclusion of 
subleading N 3 LO 3N contributions, see Refs. [TMTS] for 
first calculations along this line. In fact, such studies are 
a central goal of the recently formed Low Energy Nuclear 
Physics International Collaboration (LENPIC) [15] . 

In particular, one may expect that the unresolved 
discrepancies for spin-dependent nucleon-deuteron (Nd) 
scattering observables between calculations based on 
high-precision phenomenological NN and 3NF models 
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and experimental data at intermediate and higher en- 
ergies cannot be probed at the N 2 LO level of accuracy. 
Furthermore, recent advances in nuclear structure theory 
make it possible to extend the range of ab initio calcu¬ 
lations to mass numbers of A = 132 [2D]. The inclusion 
of 3NF contributions beyond the leading ones in few- 
and many-body calculations is the key for increasing the 
accuracy of predictions for nuclear observables and for 
systematic investigations of the role of chiral symmetry 
in nuclei and nuclear matter. We also note that correc¬ 
tions to the 3NF beyond N 2 LO (see Fig. |T]) contain such 
contributions, which have already been found earlier in 
phenomenological interaction models to be important for 
the accurate description of properties of light nuclei m- 

Numerous N 2 LO calculations of Nd scattering at low 
energy have revealed, generally, a good agreement be¬ 
tween theory and experimental data, see 0122] and ref¬ 
erences therein. On the other hand, the well-known dis¬ 
crepancies such as, e.g., the A y puzzle and the cross 
section in the symmetric space star configuration of the 
deuteron break up process could not be resolved at this 
chiral order [I] . Promising results based on chiral 3NFs at 
N 2 LO were reported by various groups in nuclear struc¬ 
ture calculations showing, in particular, sensitivity to the 
individual 3NF contributions, see Refs. is na for review 
articles. Also recent lattice simulations and no-core shell 
model calculations of light nuclei within chiral effective 
field theory demonstrate clearly the important role of the 
N 2 LO 3 NFs (see e.g., Refs. (21H28) ). Furthermore, for 
neutron-rich systems, contributions from 3NFs have been 
shown to be of critical importance for the shell structure 
and stability close to the neutron drip line (see, e.g., Ph 
1341). while for nuclear matter, saturation was demon¬ 
strated to be driven by 3NFs [35; 3D . 

While the formal expressions for the individual contri¬ 
butions of all the topologies shown in Fig. |T| have already 
been worked out 00, their decomposition in a suitable 
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FIG. 1. (color online) Different topologies that contribute to the chiral 3NF up to N 3 LO (and N 4 LO). Nucleons and pions 
are represented by solid and dashed lines, respectively. The shaded vertices denote the amplitudes of the corresponding 
interaction. Specifically, the individual diagrams are: (a) 27r exchange, (b) l7r-contact, (c) pure contact, (d) 27r-l7r exchange, 
(e) ring contributions, (f) 27r-contact and (g) relativistic corrections. See main text for details. 
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form for few- and many-body frameworks represents a 
highly nontrivial task [571439] . Due to the huge amount 
of computational resources needed for this decomposi¬ 
tion, matrix elements have been so far available only in 
a limited model space m- As a consequence, consistent 
N 3 LO three-body scattering calculations were limited to 
low energies and no studies of heavier nuclei were pos¬ 
sible. In this paper we present a novel framework that 
allows one to decompose 3N interactions in a plane-wave 
partial wave basis in a computationally much more ef¬ 
ficient way than the framework of Refs. [38) '.39] ■ This 
new method makes explicit use of the fact that all (un¬ 
regularized) contributions to chiral 3NFs are either local, 
i.e. they depend only on momentum transfers, or they 
contain only polynomial non-local terms. 

In Section m we derive the new framework for decom¬ 
posing local 3NFs efficiently in a momentum-space par¬ 
tial wave basis. In Section |III| we apply the calculated 
matrix elements of chiral 3NFs up to N 3 LO to nuclear 
matter and 3 H, discuss the partial wave convergence and 
investigate the importance of the individual topologies at 
different orders in the chiral expansion. In Section IV we 


summarize and given an outlook of future applications. 


II. PARTIAL WAVE DECOMPOSITION OF 
LOCAL THREE-NUCLEON FORCES 

A general translationally invariant 3NF can be ex¬ 
pressed as a function of the Jacobi momenta p = kl ~ k2 
and q = | [k 3 — |(ki + k 2 )], where k, denote the single 
nucleon momenta (in the following equations we will first 
suppress spin and isospin degrees of freedom): 

V123 = Vi 23 (p,q,p',q')- ( 1 ) 

Here and in the following p and q (p' and q') denote 
the Jacobi momenta of the initial (final) state. For local 
interactions, however, the momentum dependence fur¬ 
ther simplifies as such forces only depend on momentum 
transfers, i.e. on differences of Jacobi momenta: 

V123 = ^123 (p' - P, q' - q) = W^(p, q). ( 2 ) 

Note that this statement refers to unregularized forces. 
Below we will apply non-local regulators to the partial- 
wave decomposed matrix elements. The regularization 
will be discussed in more detail in Section IIIII 


Generally, the decomposition of 3NFs in plane-wave 
partial waves involves the evaluation of projection inte¬ 
grals of the form 


= J dp'dq'dpdq 

xYZ, mL , (p')V; m| , (q ')Y LmL (p )Y lmi (q)^ 3 c (p, q) ( 3 ) 


for fixed values of p = |p|,g = |q|,p' = |p'|,g' = |q'| 
and the angular momentum quantum numbers. By using 
symmetries, it is possible to reduce the dimensionality of 
the angular integrals from 8 to 5. Traditional methods 
are based on a direct discretization and numerical evalu¬ 
ation of these angular integrals [23 [35]. Due to the large 
number of external quantum numbers and momentum 
mesh points such algorithms require very large computa¬ 
tional resources for calculating all matrix elements nec¬ 
essary for many-body studies. As a result, the number 
of matrix elements of chiral N 3 LO interactions were so 
far insufficient for studies of nuclei and matter. However, 
it is possible to evaluate the basic function F defined in 
Eq. <§ in a much more efficient way by explicitly mak¬ 
ing use of the local nature of the 3NFs. Indeed, using 
rotation invariance of the potential V /23 we can write it 
as a function of three independent variables: 


El23 (p, q) = V123 (P, ?> cos 0pq), ( 4 ) 


where 


COS 0pq — 


p q 

pq 


p=|p|, 9=|q|- 


( 5 ) 


This already shows that the original eight-dimensional 
integral contains actually only three non-trivial integra¬ 
tions. The other five integrations, after employing some 
integral transformations, which are described in the ap¬ 
pendix, can be performed analytically. The final result 
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is given by 

F^r L ' mi '(p,q,p', q ') 

— X f'l'l m L +m v 

u mL—m L / ,ra z / —mi\ ± ) 

nin(Z/ +L,l' +Z) s 7 T—m L /+m L ^,1— m^+mi 
\ > L'—m L /L itilI'— mi/lmi 


,„ 2 ( 2 tt ) 4 

pp’qq’ 


\ E 

— L\,\l'—l\) 
rp'+p rq+q 

x / dpp / dqq 

J\p'-p\ J\q'-q\ 


2Z + 1 


X y l L>L(P e z +P,p) 

X yjn{qe z + q, q) 

x j d cos 6> P q/-t(cos 0 pq ) (<?, p, cos d pq ) (6) 
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with 


tfjf(a,b)= £ C^ a4mt y™“(a)l^(b). (7) 


TABLE I. Dimension and file size of the individual 3NF 
matrix element files up to N 3 LO for the different three- 
body partial waves. All matrix elements are calculated and 
stored in such a way that values of the low-energy couplings 
ci , C3 , C4 , cd, ce, Cs and Ct can be chosen freely for the dif¬ 
ferent topologies, leading to in total 12 files for each partial 
wave (see main text). For all partial waves N p = N q = 15 has 
been used. N a denotes the number of partial-wave channels 
defined in Eq. ©. All given values apply to both three-body 
parities. 


J T 

J max 

N a 

filesize [GB] 

1/2 1/2 

8 

66 

0.8 

3/2 1/2 

8 

126 

3.0 

5/2 1/2 

8 

178 

6.0 

7/2 1/2 

7 

190 

6.8 

9/2 1/2 

6 

178 

6.0 

1/2 3/2 

8 

34 

0.2 

3/2 3/2 

8 

65 

0.8 

5/2 3/2 

8 

92 

1.6 

7/2 3/2 

7 

91 

1.6 

9/2 3/2 

6 

94 

1.7 


Realistic nuclear forces also depend on the spin and 
isospin quantum numbers of the nucleons. As a com¬ 
plete basis, we choose the standard Jj-coupled three- 
body plane-wave basis of the form [3D] 

I pqa) = | pq;[{LS)J(ls)j}J{Tt)T) , (8) 


of spin-momentum operators in the form 



where L, S , J and T denote the relative orbital angular 
momentum, spin, total angular momentum and isospin of 
particles 1 and 2 with relative momentum p. The quan¬ 
tum numbers l, s = 1/2, j and t = 1/2 label the orbital 
angular momentum, spin, total angular momentum and 
isospin of particle 3 relative to the center-of-mass of the 
pair with relative momentum p. The quantum numbers 
J and T define the total three-body angular momen¬ 
tum and isospin (for details see Ref. @D]). In Eq. ^ we 
have used the fact that the 3NFs do not depend on the 
projections mj and nip , hence we omit these quantum 
numbers in the basis states here and in the following. The 
detailed basis sizes for the different three-body channels 
for our calculations of N 3 LO matrix elements are pre¬ 
sented in Table|l] Notice that as will be shown below, it 
is sufficient for our purposes to truncate the value of the 
total three-body angular momentum at J = 9/2. We 
have also verified that further increasing the values of 
•/max leads to negligibly small effects for all performed 
calculations. However, for future studies of heavier nu¬ 
clei we could extend the calculations to even larger J 
and (/ max , if necessary. 

Due to the explicit momentum dependence of the spin- 
momentum operators, the knowledge of the function F 
in Eq. © is, in general, not sufficient and the framework 
needs to be extended. This can be done in a straightfor¬ 
ward way by factorizing out the momentum dependence 


and combining the additional spherical harmonic func¬ 
tion with the ones in Eq. © by using 


y) m (x)Yi A 1 (x) = 


Z+l 

E 

L=\l-1\ 


3 2Z + 1 
4tt 2L + 1 


CmoCti ( 10 ) 


where x represents a Jacobi momentum. This strat¬ 
egy makes it possible to reduce the expressions for ar¬ 
bitrary spin-dependent interactions to the expression for 
spin-independent interactions times some momentum- 
independent spin operators. This step has to be per¬ 
formed for each momentum vector in the spin-momentum 
operators. Obviously, the efficiency of the present algo¬ 
rithm decreases with each additional su m o ver the quan¬ 
tum numbers p and L in Eqs. © and ©. Note, how¬ 
ever, that each of these sums contains only three terms 
at most. 

In order to factorize the momentum, spin and isospin 
space, for practical calculations we first calculate the in¬ 
teraction matrix elements in a L^-coupled basis: 


I PQP) = I pq; [(Ll)C(Ss)S} J(Tt)T) 


( 11 ) 


and recouple only at the end to the Jj-coupled basis 
defined in Eq. ©. Each time the factorization in Eq. © 
is applied, the spin matrix element effectively becomes 
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dependent on the quantum number /z, i.e. generally the 
matrix element in spin space can be written in the form 


(Ss)Sm s \d cr ({fj,i})\{S's')S l ms') , 


( 12 ) 


where the index i counts the number of momentum vec¬ 
tors in the spin operator. In the same way the function 
F in Eq. ^ becomes a function of the quantum num¬ 
bers Hi, i.e. it takes the form ^g 

explicit, if we consider, e.g., the case x = p in Eq. & 
the function F takes the form 


rptn L mim L / m l / fj, 

r LIL’V 


,fl (p,q,p',q') =p 


L+l 

E 

L=\L-1\ 


2L + 1 
2L + 1 


yfiO rtLm L +M j^rriLmi m L r m l i / / /x /,ol 

Xi 'LOW L Lm L lii *LlL'V G-H 

where we included the factor \J~^-p from the spin oper¬ 
ator factorization in Eq. © in this function. 

For an efficient calculation it is important to note that 
all quantities that depend on the projection quantum 
numbers m and /z are momentum independent. Hence, it 
is advantagous to factorize this dependence in the func¬ 
tion F. Specifically, for the example shown in Eq. (13) 
we can write 


form of the interaction. For example, the matrix ele¬ 
ments of the chiral long-range interactions at N 2 LO pro¬ 
portional to the couplings c\ and C 3 can be calculated 
with speedup factors of greater than 1000. Practically, 
that means that it is now possible to calculate the ma¬ 
trix elements of all interaction terms listed in Table 0 up 
to N 3 LO on a local computer cluster for sufficiently large 
basis sizes for studies of few-nucleon scattering processes, 
light and medium mass nuclei as well as nuclear matter. 
Obviously, the efficiency of the present method decreases 
with each additional internal sum in Eq. |9| and Eq. (10 1 . 


However, for all 3NF topologies up to bPLO except the 
relativistic corrections, we achieve speedup factors of typ¬ 
ically > 100 . 

Even though the present algorithm makes explicit use 
of the local nature of the chiral 3NFs, it is also possible to 
treat polynomial non-local terms. This is of immediate 
practical importance since the relativistic corrections at 
N 3 LO have precisely this form ®j. Consider, for example, 
a non-local momentum structure in the center of mass 
frame of the type 

(k 3 + k 3 ) • (k 3 - k' 3 ) = (q' + q) • (q' - q). ( 16 ) 

Such terms can be treated by factorizing the momentum 
dependence like in Eq. (|9|, for example: 


r LIL’V <Z) y ) q I — u m L -m L i ,m,/ -mi\ 1 ) 

+mi 
ilmi 


+mi 

/ j ^L' — m L /L iyllI'— m t rl 


l 


L 


(14) 


For general interactions, the function F depends on mul¬ 
tiple quantum numbers Li, hence the function takes for¬ 
mally the form F l ^)\, CPj < 7 ,lA q')- Using this decomposi¬ 
tion we can first precalculate all sums over the projection 
quantum numbers m and jn an d prestore the result in a 
function of the form Then the final matrix el¬ 

ement in LS'-coupling can be calculated very efficiently 
via: 

(pq/3\Vi2 3 \p'q'l3') = EE Flilh 1 (Pj p' j q ')> 

i {£,} 

(15) 

where values of the quantum numbers L, L', l and V are 
speficied by the PS'-coupling partial wave indices /3 and 
/?' (see Eq. (|TT])). 

Note that by deriving Eq. ( fl5| the original problem of 
calculating numerically a 5-dimensional integral for each 
matrix element as in Eq. ([3| has been reduced to the 
evaluation of a few discrete sums. The calculation and 
prestorage of the matrix elements of F l ^} v {p,q,p',q') 
can be performed relatively efficiently since only three 
internal integrals have to be performed numerically. The 
exact speedup factor of the present method compared to 
the conventional approach (55; ;35] depends on the num¬ 
ber of internal sums over Hi and Li, i.e. on the specific 


q ■ q' = qq' Y E y in (q) y *M 2 (q') e n ■ e ^ (17) 

M1 jM2=—1 

and then following exactly the steps like after Eq. ©■ 
Obviously, the algorithm becomes less efficient for non¬ 
local interactions, but the current framework turns out 
to be still more efficient than the conventional approach 
for the relativistic corrections to chiral 3NF at N 3 LO. 

Generally, 3N interactions can be decomposed in terms 
of Faddeev components: 


^123 = 5>123, 
2=1 


(18) 


whereas each of the three Faddeev components U 123 is 
symmetric in the particle labels j,k ^ i £ {1,2,3} and 
the components are related via permutation transforma¬ 
tions: 


T/( 2 ) _ p T/W p— 1 

r 123 — -*123 r 123-*123; 


= PinV&Pul. (19) 


HU p— 1 


Here P 12 3 (P 132 ) are the permutation operators that 
permute three particles cyclically (anti-cyclically) (see 
Ref. [ID])- The decomposition (18) does not uniquely 
define the Faddeev components, and the specific val¬ 
ues of matrix elements are generally convention- 

dependent. However, if evaluated between antisym¬ 
metrized wave functions, all choices and all Faddeev com¬ 
ponents give the same results. In contrast, the matrix 
elements of the completely antisymmetrized interaction, 


(1 + P 123 + P 132 ) U 123 (1 + -P 123 + -^ 132 )) 


( 20 ) 
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are unique. Since the application of the permutation op¬ 
erators in Eq. (20) in a momentum partial wave basis is 
non-trivial and can induce numerical uncertainties (see 
Ref. [35J), we calculate and provide matrix elements of 


both (pqa\Vi 23 \p'q'cx') and (pqa\Vi 2 3 \p'q'a'). This of¬ 
fers the flexibility to transform the matrix elements of the 
Faddeev components to a harmonic oscillator basis and 
perform the antisymmetrization directly in this basis. 
The resulting matrix elements can then be used for cal¬ 
culations within e.g. the no-core shell model, the valence- 
shell model, the coupled cluster or the self-consistent 
Green’s function framework. On the other hand, the 
antisymmetrized matrix elements (pgo'|V r 1 ^ s 3 |p , g , a ’) can 
be directly used for studies of infinite nuclear matter 
or few-body scattering processes. Due to the large ba¬ 
sis sizes shown in Table ]T| the uncertainties induced by 
the antisymmetrization are very small. We have checked 
that for three-nucleon systems the obtained binding ener¬ 
gies based on the Faddeev components and the antisym¬ 
metrized interactions are identical within the sub-keV 
level. 

For the application of the permutation operator P 123 
in the three-body plane-wave basis defined in Eq. ([ 8 ]) it is 
important to implement this operator in an efficient and 
numerically stable way. This is in particular important 
for calculations in large bases such as those shown in Ta¬ 
ble [TJ Several different expressions have been derived for 
the permutation operator (see, e.g., [301 [33] ). All these 
expressions suffer from problems due to ratios of possi¬ 
bly very small numerators and denominators, which can, 
lead to numerical instabilities for partial waves with large 
angular momenta. For our calculations, we use a novel 
improved implementation. Following the derivations of 
Ref. [30], it is straightforward to derive the following ex¬ 
pression: 


{pqa\ P123 \p'q'oi') = E \!J 3 J' 3 '& 
c,s 


f L S J W V S' J' 1 

x { l 1/2 j } < V 1/2 f \ 

( C S J j [ C S J j 


1 S \ 

V 1 \ 1/2 S S'j 

x(-l fVff 1 * { 1,/2 1,/2 T \ 

' [1/2 r t' j 


x8tt / dcos8. 


$(p' - |l/2p + 3/4q|) S(q' - |p - l/2q|) 


pq 


n'2 


P* Q 

x E (P(—l/2p^3/4q, p—T/2q). (21) 


The key difference to other expressions is the fact that 
we directly perform the angular integrals as shown in 
Eq. (21) without decomposing the angular dependence 
of the spherical harmonic functions any further. This 
implementation turns out to be numerically more efficient 


and, most importantly, is perfectly stable even for large 
values of angular momenta. The matrix elements of the 
operator P 132 are also given by Eq. plj ) (see Ref. [30] 1. 

As a first application, we use the new framework to 
calculate the matrix elements of the chiral 3NF up to 
N 3 LO. The individual topologies are shown in Fig. [TJ In 
addition to the pion decay coupling F n and nucleon axial- 
vector coupling constant gA , the chiral 3NFs up to N 3 LO 
depend on 7 LECs in total, namely Ci, C 3 , C 4 , cd, ce, Cs 
and Ct [ZHS 02] , for which we want to keep the flexibil¬ 
ity of being able to change their values after performing 
the partial-wave decomposition. Specifically, the indi¬ 
vidual topologies ar^j](see also Fig. I]): (a) 27 t exchange 
[ci,C 3 ,C 4 ,1], (b) l7r-contact [cp], (c) pure contact [ce], 
(d) 27T-17T exchange [ 1 ], (e) ring contributions [ 1 ], (f) 2tt- 
contact [Ct] and (g) relativistic corrections [Cs,Ct,1]- 
Here, the couplings in square brackets denote those free 
LECs which appear in each of these topologies, whereas 
the constant ’F indicates contributions that do not con¬ 
tain any of the indicated 7 LECs. This leads in total to 
12 individual contributions for each three-body partial 
wave as defined in Eq. 

For each of these contributions we first compute one of 
the Faddeev components {pqo\V \23 \p'q'a') for all partial 
wave channels and basis sizes shown in Table |T] Subse¬ 
quently, for the calculation of the antisymmetrized ma¬ 
trix elements (pqa | | p'q'a') we apply the permutation 

operator as defined in Eq. ( |2l| ). The application of the 
permutation operator in the partial wave basis involves, 
in principle, a complete sum over intermediate partial 
wave quantum numbers (see, e.g., Ref. [39]). Thanks to 
the large basis sizes shown in Table [T] we ensure that 
the antisymmetrization leads to well-converged results. 
Furthermore, it is straightforward to generate arbitrary 
other products of the permutation operator P 123 and the 
Faddeev components in a computationally very efficient 
way. Such products appear, for example, in Faddeev 
equations for few-body scattering problems m- 


III. APPLICATION OF CHIRAL 3N FORCES 
AT N 3 LO TO NUCLEAR MATTER AND 3 H 

Figs. [2] and [3] illustrate the convergence of the par¬ 
tial wave expansion of the calculated chiral 3NFs. These 
figures show the 3NF contributions to the energy per 
particle of neutron matter (Fig. 2) and symmetric nu¬ 
clear matter (Fig. 3) in Hartree-Fock approximation for 
the individual partial-wave channels and 3NF topologies. 
These results provide direct insight in the required num¬ 
ber of partial-wave channels in mean-field calculations. 
Of course, the number of required partial-wave channels 


1 We correct a misprint in Eq. (4.14) of Ref. [9] and use the param¬ 

eter values /3 q — / and = /, corresponding to the ” minimal 

nonlocality” choice of the potential (see Ref. m- 
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FIG. 2. (color online) Partial wave contributions to the energy per particle to neutron matter in the Hartree-Fock approximation 
at nuclear saturation density (kp = 1.7fm _1 ) for the individual 3NF topologies up to N 3 LO. For the shown energies we use 
the coupling constants Cs = Ct = 1 and d = lGeV -1 . All results show the accumulated energy contributions, including 
all contributions up to the given partial-wave channel. Here P denotes the three-body parity P = (— l) L+l and J is the 
three-body total angular momentum, as defined in Eq. J8|. In neutron matter only channels with T = 3/2 contribute. The 
exact benchmark results are calculated following Refs. 11711181 . 
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FIG. 3. (color online) Partial wave contributions to the energy per particle to symmetric nuclear matter in Hartree-Fock 
approximation at nuclear saturation density (kp = k p F = 1.35 fnG 1 ) for the individual 3NF topologies up to N 3 LO. For the 
shown energies we use the coupling constants Cs = Ct = 1 and a = 1 GeV -1 . The results for the 27r-contact topology are 
scaled by a factor of 1/4 for presentation purposes. All results show accumulated energies, including all contributions up to 
the given partial-wave channel including both three-body parities. J and T are the three-body total angular momentum and 
isospin, as defined in Eq. Q. The exact benchmark results are calculated following Ref. 118} . 


has to be checked in more detail in realistic nuclear struc¬ 
ture calculations, for which also many-body correlations 
are important. Moreover, these results also serve as a di¬ 
rect benchmark of the calculated matrix elements, since 
the energy contributions in the Hartree-Fock approxima¬ 
tion have been already calculated independently based 
on the operator expressions directly, see Refs. [T71[TH] for 
details. 

Specifically, for the results shown in Fig. [2] we use kp = 
1.7 fm -1 for the neutron Fermi momentum, which corre¬ 
sponds to a neutron number density of n ~ 0.166fm -3 . 


Since the Fermi momentum serves as an ultraviolet cut¬ 
off, we do not need to regularize the interactions for these 
benchmark calculations. In neutron matter only matrix 
elements in the three-body isospin channel T = 3/2 con¬ 
tribute, whereas the N 2 LO topologies that include the 
low-energy coupling constants C 4 , cp and cp vanish ex¬ 
actly [33j. The detailed expressions for the Hartree-Fock 
energies for neutron matter in terms of the antisym¬ 
metrized matrix elements V^ 3 are given in Ref. • The 
results of Fig. [2] show that matrix elements up to the 
partial wave channel with J = 5/2 and both three-body 
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FIG. 4. (color online) Contributions of the individual topologies to the triton energy for three different NN interactions of 
Refs. IIIII 22 I . The LECs cd and ce are chosen to be of natural size and are fitted to the experimental triton binding energy. The 
topologies ci, C 3 and C 4 contain contributions from N 2 LO and N 3 LO, Cn and ce are pure N 2 LO contributions and 27r(no Ci), 
27r-l7r, rings, 2-7r-cont and 1/m corrections are N 3 LO contributions. See main text for details. 


parities P = (— \) L+l can provide significant contribu¬ 
tions to the energy, whereas higher partial waves give 
only small corrections. Overall, we find that the results 
including all contributions up to J = 9/2 are very well 
converged and show excellent agreement with the exact 
results (see also 0Sj). 

For symmetric nuclear matter we find a very similar 
convergence pattern. In contrast to neutron matter here 
all 3NF topologies shown in Fig. |T] and also both three- 
body isospin channels, T = 1/2 and T = 3/2, contribute 
to the energy. For the results shown in Fig. [3] we fix 
the neutron and proton Fermi momenta to k F = k p F = 
1.35 fm _1 , which again corresponds to a total number 
density of n ~ 0.166 fm“ 3 . We show the contributions 
to the energy for the individual partial-wave channels, 
whereas here each \J , 7~] channel includes contributions 
from both three-body parity channels P = ±1. Again, 
we observe excellent partial wave convergence and es¬ 
sentially perfect agreement with the exact Hartree-Fock 
results. 

Next, in Fig. [3] we illustrate the contributions of the 
individual topologies to the binding energy of 3 H. In 
order to probe scheme dependence of our results, the 
present calculations use three different NN interactions: 
the N 3 LO potentials of Ref. [TSj with the cutoff combina¬ 
tions A/A = 450/500MeV and 550/600MeV (EGM) and 
the N 3 LO potential of Ref. [IT] (EM). For our calcula¬ 
tions we fix the values of the LECs ci, C3, C4, Cs and Ct 


consistently to their values of the NN interactions. No¬ 
tice further that the N 3 LO corrections to the two-pion 
exchange topology involve contributions which account 
for finite shifts in the LECs Cj in the N 2 LO 3NF ex¬ 
pressions [§]. All these effects are properly taken into 
account. The couplings cd and ce are fixed to the ex¬ 
perimental binding energy of 3 H and by requiring the 
values of those LECs to be of natural size. For each 
of the NN potentials we used three sets of values, their 
specific values for the different interactions are given in 
the legend of Fig. [4] The matrix elements of the 3NF 
are regularized by applying a non-local regulator of the 
form f R (p , q) = exp[-((p 2 +3/4q 2 )/A^ N ) 3 ] as in Ref. gg], 
whereas the value of the cutoff A 3 N is chosen to be consis¬ 
tent with the cutoff value A of the corresponding NN in¬ 
teraction. Note that all results of this figure are based on 
the consistent wave functions for each of these fits, which 
obviously complicates a detailed quantitative compari¬ 
son of results based on different fits or NN interactions. 
However, we do not expect that this affects the main 
qualitative features of the results, which can be sum¬ 
marized as follows: first, the size of the contributions 
of the individual topologies depends sensitively on the 
employed NN interaction as well as on the values of the 
LECs cd and ce■ Second, the contributions of topologies 
at N 3 LO are not suppressed compared to those at N 2 LO 
for the presently employed 3N interactions. Similar find¬ 
ings have been found before for 3 H using a restricted 


















number of partial waves m and for the mean-field con¬ 
tributions of the N 3 LO 3NFs to neutron matter and sym¬ 
metric matter mm- Third, the perturbativeness of the 
3NFs also depends sensitively on the employed NN inter¬ 
action. While, e.g., for the interaction EGM 550/600 the 
wave function is not strongly affected by variations of the 
short-range couplings, for the other two NN interactions 
the contributions from the long-range 3N topologies also 
change when the LECs cn and ce are varied. 

When interpreting these results, it is important to keep 
in mind that neither the individual nor the total contri¬ 
bution of the 3NFs to the binding energies of nuclei are 
experimentally observable. The considered contributions 
of individual topologies represent, strictly speaking, bare 
quantities while all estimations based on power counting 
refer to renormalized ones, see also Ref. m- Moreover, 
the employed nonlocal regulators do not actually cut off 
all short-range pieces of the 3NFs, so even the long-range 
3NF topologies do contain short-range contributions af¬ 
ter regularization (see Ref. m for a related discussion). 

Clearly, expectation values of such short-range ad¬ 
mixtures are strongly scheme dependent. Recently, a 
novel type of chiral NN interactions has been devel¬ 
oped [m uni an. For these new interactions, the long- 
range topologies are regularized by local regulators that 
act only on the relative particle distance in coordinate 
space and allow for a cleaner separation of long- and 
short-distance physics. Indeed, in Ref. [50] it was shown 
explicitly that 27r-exchange contributions dominate at 
long- and intermediate distances when local regulators 
are employed for NN interactions. 


IV. SUMMARY AND OUTLOOK 

In summary, the new framework presented in this work 
makes it possible to include the chiral 3NF at N 3 LO and 
beyond in ab initio few- and many-body calculations. For 
non-local regulators, the calculated matrix elements can 
be immediately used within various many-body frame¬ 
works for novel studies including all NN and 3N contri¬ 
butions consistently beyond N 2 LO. The application of 
the regularization scheme employed in Ref. [TOl US] for 
the NN forces to 3NF will require a generalization of the 
present framework. In contrast to the considered non¬ 
local regulators, which only depend on the magnitude of 
the Jacobi momenta and act just as multiplicative fac¬ 
tors in momentum space, local regulators do mix differ¬ 
ent partial waves. The application of these regulators in 
momentum space can, e.g., be implemented using their 
partial wave decomposed form and performing numerical 
evaluation of the corresponding folding integrals. Work 
along this line is in progress. 

We also emphasize that a careful investigation of sys¬ 
tematic uncertainties in few- and many-nucleon calcu¬ 
lations represents an important goal in nuclear physics, 
see also discussions in Refs. hsied. The availability of 
matrix elements of NN and 3N interactions at different 


orders in the chiral expansion and within different regu¬ 
larization schemes will provide an important contribution 
towards a better understanding of the systematic uncer¬ 
tainties and ultimately allow for precision tests of chi¬ 
ral dynamics in nuclear systems. Last but not least, we 
emphasize that the presented framework can be straight¬ 
forwardly applied to 3NF at order N 4 LO [52] [55] and to 
EFT interactions with explicit A-degrees of freedom j5D . 
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Appendix A: Appendix: 

Angular integrations in partial wave decomposition 

In this Appendix we describe the integral transforma¬ 
tions which allow for a decoupling of the three non-trivial 
integrations over spherical harmonics in the partial wave 
decomposition of the 3NF from the other five which can 
be performed analytically. Let us start with the Eq. ([3]) 
and add a radial integration over p' and q' in order to 
have a translationally invariant measure. We can achieve 
this by introducing additional integrations with delta- 
functions via 

F ?ii7 mL ' mi ' (P’V’P' ’J) = J d 3 q'd 3 p'dqdp 

x6(p’ - p 1 )6(q’ - q 1 )Y£, mLl (£' )Y l J mv (4' )Y LmL (p )Y lmi (q) 

xVmlp' - p,q' - q)- (Ai) 

Now we can make a substitution 

p' -> p' + p and q' A q' + q, (A2) 
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which leads to 


F ™l™r L,rni1 {p,q,p',q') = J dVdVdqdp 

x5(p' - |p' + p|)%'- |q' + q|)V 1 1 g(<f,p',cos0~,q) 
xY L'm L , (p' + P )Y l r mi , (q' + q )Y L m L (p )Yi mi (q). (A3) 


In this form it is manifest that the integrals factorize and 
we need to calculate 

p/ 2 q/2 J dqdpdp' J d<j)q’8{p' — |p' + p| )8{q' — |q' + q|) 

xY l, mLl (p' + p )Y l 1 mv (q' + q )Y LmL (p )Y lmi (q), (A4) 

where 4>q' is an azimutal angle of q'. This integral does 
not involve the 3NF and can be calculated analytically. 
The full result is given by Eq. 
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